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A literature review is presented of the work done so far in the numerical modeling of 
wind flows. Pertinent computational techniques are described, as well as the necessary 
assumptions used to simplify the governing equations. A steady-state model is outlined, 
based on the data obtained at the Deep Space Communications complex at Goldstone, 
California. 


I. Introduction 

In recent years, due to the increasing cost of oil and other 
fossil fuels, wind-generated electricity has been actively sought 
as a possible alternate and renewable energy source. At the 
present time, there exists a national R&D program to test the 
engineering and economic feasibility of wind-driven turbines. 

In determining the site and construction of such turbines, a 
necessary starting point is the collection of data on wind 
speed, direction and duration. Since this data collection is 
hmited in scope due to resource restrictions, alternative means 
must be sought. This siting problem can be stated as follows: 
for a certain area, given the topography of the terrain and 
given the wind direction and speed at one or more of the area 
boundaries, it is desired to predict the wind distribution at all 
points within the stated boundaries. 

The solution of the problem can be quite complex due to 
three main sources of difficulty. First, the problem is 
complicated by the topography of the terrain, which may be 
highly irregular, requiring a large amount of data to describe it. 
Second, the wind data, given at one or more of the boundaries, 
is highly irregular with respect to temporal variation in 
magnitude and with changes of direction. Therefore, there are 


difficulties related to processing the large amounts of data 
necessary for describing these boundaries. Third, and of equal 
importance, is the fact that the fluid dynamics and energy 
equations describing the wind flow are too comphcated to be 
solved by analytic means in closed form. Because of these 
reasons, any analysis for the prediction of wind velocities over 
an arbitrary terrain must be treated by numerical, i.e., 
computer methods. 

This article is written with the aim of introducing the 
reader to various aspects of wind prediction research. The first 
section presents the results of a literature search made during 
June — December 1979 into this subject and it also describes 
the past DSN activities in this field. 

The next four sections discuss the physical and computa- 
tional aspects of pertinent numerical modeling. The article 
concludes with suggestions for future work. 

II. Literature Review 

During the period from June — December 1979, a review of 
the literature pertaining to wind velocity predictions has been 
conducted in the Energy Group of the DSN Engineering 
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Section. It became apparent that most of the research effort in 
this area had been concentrated at only a few laboratories and 
universities. The results obtained at these centers are presented 
below. 

A. Work Performed at Los Alamos Scientific 
Laboratory 

The Los Alamos Scientific Laboratory (LASL) has a strong 
tradition of computer modeling dating back to the 1940s. The 
laboratory possesses some of the most modern computers 
available and an excellent technical staff. This combination has 
produced several classic techniques for fluid flow modeling 
such as the Marker And Cell (MAC) (Ref. 1) and the Particle 
In Cell (PIC) methods (Ref. 2). In recent years, LASL has 
published a number of papers dealing with computer programs 
that model the spread of wind-borne pollutants (Refs. 3-5). 
These programs are characterized by (1) solutions of the 
complete, unsteady Navier-Stokes equations, and (2) the appli- 
cation of these solutions to flows past obstacles with a regular 
geometry, e.g., cubes and rectangular blocks. The results 
obtained show very good agreement with experimental obser- 
vations of pollutant dispersion. 

A penalty paid for obtaining such good results is the high 
cost incurred by operating these programs. Nearly all the 
storage of the 64,000-word fast core memory of a CDC — 
7600 computer is taken up by these codes. Within these 
programs, there are approximately 3300 computational cells, 
each of them requiring approximately 1 .5 sec of computation 
time per cycle; this means significantly more than 1 hour of 
computation time for each cycle. Some of the simpler 
programs developed at LASL (such as SOLA) still require 
large amounts of computer time and cost. A rough estimate 
(Ref. 10) shows that these costs are proportional to (N^), 
where N is the number of mesh cells into which the length 
dimension is divided, and with (/f^), where K is proportional 
to the distance that the mesh boundaries must be kept away 
from the flow region of interest. Therefore, computational 
costs can be kept within some reasonable limits by accepting 
the minimum resolution that will still yield accurate results, 
that is reducing A^ and providing realistic boundary conditions, 
which means minimizing /f. It was shown (Ref. 6) that these 
cost penalties still remain high even when the physical 
assumptions of the problem are slightly relaxed. 

B. Work Performed at Lawrence Livermore 
Laboratory 

A second center which is involved in research on wind field 
prediction is the Lawrence Livermore Laboratory (LLL). The 
impetus for their studies is also the desire to predict air 
pollution dispersion. The papers published by LLL point out 
both the diversity which can exist in computer modeling and 


the marginal results that one often encounters when such 
numerical modeling is applied to wind prediction. In contrast 
with the LASL efforts, the model used at LLL (Ref. 7) 
employs variational techniques to obtain a time-independent, 
three-dimensional wind field model. Despite the sophisticated 
techniques involved, the comparison of numerical predictions 
with terrain measurements is not good. “Typically, 60% of the 
time the calculations of pollutant dispersion were within an 
order of magnitude (of terrain measurements).” The chief 
sources of error were found to be, in decreasing order of 
magnitude, (l)wind direction, (2) topography, (3) diffusion 
parameters, (4) source strength, and (5) wind speed. Probably 
because of its time-independent characteristics, the LLL 
program is shorter to run than corresponding LASL models; 
the LLL program takes 80% of the large core memory on a 
CDC — 7600 computer and uses 2 to 3 minutes of CPU time 
to generate approximately 30,000 grid points above the 
terrain. There is no documentation as to how many iterations 
(i.e., sweeps) are needed to achieve the steady-state conditions. 

C. Work Performed at Colorado State University 

Two university centers have also done important work 
regarding different aspects of wind modeling. Research done at 
the Colorado State University, Fort Collins (CSUFC), has 
centered on experimental wind tunnel modeling of atmo- 
spheric flows over various topographies. Specialized meteoro- 
logical wind tunnels are used in these studies (Refs. 8 and 9). 
The wind tunnel at Ft. Collins is capable of simulating 
thermally stratified atmospheric boundary layers. This type of 
modeling is also being used in France, Australia, and New 
Zealand (Refs. 10-12). The foremost advantage to this 
approach is that it lowers costs if the terrain conditions can be 
reproduced accurately in a wind tunnel. Also, it becomes 
much cheaper to obtain pertinent data from wind tunnel 
measurements rather than from instrumentation erected and 
monitored on the actual terrain. The results published by 
CSUFC indicate that there exists enough similarity between 
wind tunnel patterns and actual wind fields to justify this 
approach (Ref. 13). It is still unknown how adequate this 
modeling can be in non specialized or nonmeteorological wind 
tunnels, and very little information about it has so far 
appeared in the literature (Ref. 14). The CSUFC interest in 
computer modeling seems to be centered on the analysis of 
atmospheric turbulence. Significantly, the results which were 
obtained from numerical models were then checked with wind 
tunnel measurements (Ref. 60). 

D. Work at Other Universities and Centers 

The other university center which has done long-range, 
sustained research on wind fields is Pennsylvania State 
University (PSU), where H. A. Panofsky and his co-workers 
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have conducted a number of investigations on the effect of 
terrain roughness on wind profiles (Refs. 15-19). Their data is 
probably the best record on how different types of terrain 
roughness such as forests, cropland, and desert can affect the 
wind profiles. How important this information would be in a 
numerical model of wind prediction, is not known at the 
present time. Work has also been done at Sandia Laboratories, 
Albuquerque, New Mexico on numerical wind modeling (Ref. 
58), but the present status of this research is unknown. 


The NASA-Lewis Research Center has been involved in 
designing and building their first lOO-kW(e) wind propeller and 
working in cooperation with the Department of Energy as 
consultant in implementing similar or larger units at various 
other locations. However, relative to numerical modeling of 
wind velocities, it was not evident that NASA-Lewis was 
involved in detail. 


E. JPL Activities 

The wind prediction studies at JPL were aimed at determin- 
ing the feasibility of using wind power as part of an energy 
system to supply the energy needs for the Deep Space 
Network at Goldstone, JPL in Pasadena, and Edwards sites. In 
1974, a Wind Power Feasibility Study committee at JPL 
published a preliminary report (Ref. 59), which contained a 
section related to wind siting. The report limited itself to 
discussing observed wind speeds at different Cahfornia sites. 
This data was based on a comprehensive report of observed 
wind speeds at 137 California sites (Ref. 42). In parallel to 
that effort, another program was undertaken to determine the 
effects of wind loading on the structure of the tracking 
antennas at Goldstone. For the latter program, continuous 
records of wind speed and direction were obtained at six 
Goldstone sites (including different heights at the same site) 
during the period from October 1974 through July 1976. 
These records form one of the most comprehensive wind logs 
available, despite considerable gaps, due mainly to equipment 
malfunctions. Besides the terrain measurements, a private 
consulting firm (MRI)‘ developed a computer model for wind 
prediction at the Goldstone site. The results of this program 
(Ref. 21) show it to be quite unreliable and very expensive at 
about $300/run. A possible cause for this poor performance is 
given as due to “lack of sufficient fineness in the rectangular 
grid representation of the terrain.” At the present time, all of 
the data taken by the Goldstone wind energy measurement 
system has been processed and a report is expected to be 
released shortly. 


F. Conclusions of Literature Search 

From the literature search briefly described above, one can 
draw the following conclusions: 

(1) No numerical model presently available can predict 
with sufficient accuracy wind speeds and directions 
over an arbitrary terrain. 

(2) Building such a model is feasible; this has been shown 
by the computer models developed at LASL for 
pollutant dispersion. 

(3) The results of any such model can be checked cheaply 
and accurately with the aid of wind tunnel experi- 
ments. 

(4) JPL has one of the best data bases available to serve as 
an input for such a model. This is extremely important 
since a computer program is only as good as the inputs 
that the program starts with. 

Building such a numerical model in a simple way would be of 
great help to the DSN in the decision of wind turbine locations 
to achieve a sizable reduction in the energy consumption and 
cost. The sections that follow address themselves to the 
“mechanics” that enter into writing such a new program. The 
fluid mechanics, mathematics, and computational bases for 
such a model are discussed in a tutorial manner starting from 
elementary notions and briefly touching some state-of-the-art 
techniques. 


III. Mathematical and Physical Description 
of the Probiem 

The flow dynamics of many liquids and gases can be 
described by means of a mathematical expression known as 
the Navier-Stokes equation. This is a vector equation which 
relates the forces and accelerations acting on a fluid particle. 

In cartesian coordinates, the x, y, and z components of this 
equation are respectively; 
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The first term on the left hand side of these equations 
represents a temporal acceleration. The remaining terms on the 
left hand side represent accelerations due to the motion of the 
fluid (convective accelerations). The right hand side of the 
equations contains the force terms due respectively to pres- 
sure, viscous and bulk effects. 


exist for simplified, linearized forms of the Navier-Stokes 
equations (Ref. 22). These simplified forms occur when the 
nature of the flow allows the investigator to neglect some of 
the terms in comparison with other terms. For example, for a 
flow which has the same axial velocity profile at all locations 
on the Z-axis, we can neglect terms that contain 9F /,8z. 

It is only in the last decade that the advent of large 
computers has made possible the solution of the full Navier- 
Stokes equations by numerical methods. In considering such 
methods, one needs to know the answers to the following 
questions: 

How are the equations approximated for numerical 
processing? 

What are the numerical errors and instabilities that can 
arise? 

What techniques are available to solve the flow problems of 
interest by fast, inexpensive, and accurate means? 

Some answers to these questions are considered next. 


IV. Numerical Procedures 


The three Navier-Stokes equations contain five unknowns: 
the velocity components V^, V^, and V^, the pressure p and 
the density p. Even if the problem is such that p can be 
considered as a known constant, there will still be four 
unknowns. In order to have as many equations as there are 
unknowns, an additional equation is needed. The needed 
relation is the conservation of mass, known as the continuity 
equation. The continuity equation for an incompressible fluid 
is: 


A. Finite Differences and Finite Elements 

There are two main numerical techniques available for the 
solution of fluid flow phenomena. They are the finite- 
difference method and the finite-element method. In the first, 
the differential equations to be solved are transformed into a 
set of difference equations. Difference equations are algebraic 
in nature. They consist of additions and subtractions of 
variables at discrete points, in place of the continuous 
derivatives that occur in differential equations. 


dV dV dV 
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Hence, for the case where p is constant (incompressible flow), 
the continuity and the three components of the momentum 
equations (Navier-Stokes) will form a well-posed system of 
four equations and four unknowns. 

An important characteristic of the Navier-Stokes equations 
is the nature of the convective acceleration terms. These 
accelerations consist of dependent variables (velocities) multi- 
plied by their derivatives, for example multipUed by 
bVJbx. Such terms are called “nonlinear”, and the equations 
that contain these types of terms are called nonlinear 
equations. Nonlinear equations are very difficult to solve by 
analytic methods. A limited, if important number of solutions 


In the second technique, the finite element method, 
variational calculus is used to solve the differential equations 
of the problem converting it to an integral form and following 
a minimization path. First, the region of interest is divided 
into discrete elements. The elements are assumed to be 
connected only at the common nodes of their boundaries, (for 
example, nodes a, b, c, in Fig. 1). The finite-element method 
searches for the values of the dependent variables at all these 
common nodes within the domain. The variational techniques 
which are used (for example, the Bubnov-Galerkin method) 
consist of setting up integrals of the differential equations and 
then minimizing the values of these integrals at each of the 
nodes in the domain. 

There is no clear-cut choice as to which of these two 
methods is better for the solutions of fluid mechanics 
problems. Finite differences are easier to set up than finite 
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elements. In addition, there is a large body of experience 
dealing with finite difference techniques as applied to numeri- 
cal fluid dynamics (Ref. 23). This experience has to do with 
questions of numerical stability and convergence and is very 
important when writing new numerical schemes. The main 
disadvantage of the finite-difference method is that it does not 
lend itself very easily to incorporating irregular physical 
boundaries. 


For the backward-difference scheme: 



V - V , 

n n-1 

X - X , 
n n-1 



( 6 ) 


Finally, for the central-difference scheme, we can take 
differences relating points on either side of n and obtain: 


On the other hand, finite element schemes, originally 
devised for structural and elasticity problems, have only 
recently been applied to fluid mechanics problems. Their great 
advantage is the ease with which the solution domain can be 
discretized. Therefore, irregular physical boundaries are easily 
described by this method. However, finite-elements methods 
suffer from difficulties in deaUng with nonhnear terms of the 
Navier-Stokes equations and with incorporating the incompres- 
sibility conditions. While new methods are being devised to 
overcome these difficulties (Refs. 24, 25), finite elements 
methods also seem to require shghtly more computing time 
than finite difference methods (Ref. 26). 


.. Vn.l-Vn-l Vn^l-Vn-l 
9a: 2M 

Because all difference formulas are only approximations to 
derivatives, they will contain some errors when compared with 
the actual values of the derivatives. The magnitude of the 
errors depends on the finite difference form used. It can be 
shown mathematically (Ref. 28) that of all three schemes 
presented above, the central-differences gives the most accur- 
ate result. 


It can be shown (Ref. 27) that finite differences are written 
for points in space and finite elements are spatial integrators of 
point formulas. However, no comparisons have been published 
regarding the efficiency and cost of these two methods when 
applied to the same problem. Mainly because of the difficulties 
still remaining with the finite element methods, we will 
concentrate our attention to the solution of the Navier-Stokes 
and continuity equations by means of finite-difference 
techniques. 


In general, velocities and pressures do not vary in one 
direction only. 'In rectangular coordinates, two- or three- 
dimensional computational cells, as those shown in Figs. 3a 
and 3b, are often used. 

For two- or three-dimensional flows, one has to decide how 
to assign the variables within each cell. For the two- 
dimensional case, the most popular schemes are illustrated in 
Figs. 4a and 4b. 


B. Derivatives and Computationai Meshes 

There are a number of ways by which derivatives can be 
approximated by finite differences. Suppose, for example, that 
one is interested in computing the changes in the x-component 
of the velocity vector with respect to the position on the 
x-axis, that is bVjbx. Consider a computational mesh, with 
equal mesh intervals A/i, in the x-direction, as shown in Fig. 2. 
Denote by (n-1), (n) and (n+1) a series of three adjacent 
points at the intersections of the grid lines. These points are on 
a grid line parallel to the x-axis. In general, the velocity 
component has different values at each of these points; 
these values are denoted respectively by and V^^.. 

The change in with respect to x, evaluated at the point n, 
can be approximated by the differences between the velocity 
values at these points. For example, for the forward- 
difference scheme: 


dV 

X 

0X 


V V 

n + l n 

X ^ - X 

n+1 n 


V - V 

n + 1 n 

Ah 


( 5 ) 


In Fig. 4a, averages at the center of the cell are used to 
determine fluxes through the cell’s edges. The staggered mesh 
arrangement shown in Fig. 4b is more advantageous for 
incompressible flows (Ref. 29). Here, the horizontal velocity 
component is centered on the right and left faces of the 
cell, while the vertical component V is centered on the 
bottom and top. The pressure and density are cell-center 
averaged. 

In addition to assigning the variables within each cell, one 
also has to “optimize” the computational grid size. An optimal 
grid size is one which yields good resolution for the fluid flow 
phenomena. It is natural to expect that the smaller the grid 
size, the better the resolution. However, a small grid size also 
implies a greater number of cells covering the domain of 
interest. With a larger number of cells, both computational 
time and cost will increase. Therefore, a good compromise is 
to use a variable mesh size having small cells in regions of rapid 
flow changes and larger cells in regions where the flow does 
not change drastically (Refs. 3 and 36). 
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In the next section, it will be shown that the selection of 
mesh size and the way the variables are assigned within each 
cell are part of the computational stability problem. Meaning- 
ful solutions of the Navier-Stokes equations can be obtained 
only when written in a certain form and when the finite 
differences are forced to take into account the direction of the 
flow velocities. 

C. Problems of Stability and Convergence 

When finite-differences are substituted for the derivatives, it 
may happen that the numerical solution behaves very differ- 
ently compared to the solution of the corresponding differ- 
ential equation. The main causes for such discrepancies are 
instability and poor convergence. 

Since any computer memory is hmited, the numbers that 
are stored and worked with in computers can be represented 
only by a finite number of digits. Numbers having more digits 
are “rounded-off’ according to algorithms designed into the 
machine. Any such rounding-off will contribute to errors in 
the final result. Therefore, it is possible that round-off errors 
will propagate through a computer program and, as a result, 
small changes in the initial data will be translated into very 
large changes in the output. Numerical schemes that suffer 
from such defects are said to be unstable. The growth of errors 
due to the presence of possible extraneous solutions to the 
difference equations, can also cause instability (Ref. 61). 


into Taylor series and analyzing the higher order terms. 
Although Hirf s method can be applied to nonlinear equations 
with variable coefficients, it is uncertain yet whether it can be 
applied to some of the common complex equations. 

It will be shown later that computational stability can be 
improved by a judicious choice of finite-difference expres- 
sions. In particular, a certain scheme of differentiation, called 
the “upwind” or “donor-cell” method, has been found to 
greatly improve the stability of the numerical solutions of the 
Navier-Stokes equations. However, before this method is 
described further, convergence must be looked at in more 
detail. 

As mentioned earlier, convergence is related to the accuracy 
of the numerical solution, and it is an indication of how close 
the numerical solution approximates the solution of the 
continuous differential equation. The laws of fluid mechanics, 
as given by Eqs. (1), (2), (3), and (4) can be written in 
different forms, to which correspond different numerical 
schemes. Some of these schemes have better convergence than 
others. In particular, when the momentum equations are 
written in what is called a “conservative form,” it has been 
argued (Refs. 29, 33 and 34), that the finite difference 
equations corresponding to this form give more accurate 
results than when the equations are in a “non-conservative” 
form. 


Another possible type of error has to do with the 
convergence of the algorithm used. In general, finite differ- 
ences are only approximations to the continuous derivatives 
that they replace, and therefore, the numerical solution will 
always be in error compared to the solution of the original 
differential equation. Mathematically, one defines a con- 
vergent finite difference scheme as one in which the finite 
differences solution approaches the continuum differential 
equation solution, as the mesh size approaches zero (Refs. 23 
and 35). 

The definitions of stability and convergence can provide 
only rough guidelines for practical numerical calculations. In 
particular, stability analysis is still in its infancy. The oldest 
method for testing the stability of a numerical scheme is due 
to von Neumann (Ref. 31) and is based on expressing the 
dependent variables in terms of Fourier series. When the 
Fourier components are substituted into the finite-difference 
equation, the decay or growth of each mode in these 
components will show whether the scheme is stable or 
unstable, respectively. While this technique works fairly well 
for linear equations with constant coefficients, it is inaccurate 
when applied to nonlinear equations such as the Navier-Stokes 
equations. Hirt (Ref. 32) has introduced another method, 
based on expanding each term of the finite difference equation 


Although there is no rigorous proof for this statement, and 
some exceptions have been found (Ref. 23), in general, it is 
advisable to write the momentum equations in the conserva- 
tive form (Ref. 35). 


D. Conservation: “Upwind-Differencing” 

An equation is said to be in “a conservative form” when it 
can be written as a sum of a time derivative and a spatial flux. 
For example, the one-dimensional momentum equation can be 
written in conservative form as 
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The corresponding “non-conservative” form of the momentum 
equation is; 
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Note that Eq. (8) can be obtained from Eq. (9) by adding to 
the latter the continuity equation multiplied by The 
corresponding finite-difference equations for Eqs. (8) and (9) 
are, respectively, 
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In the above equations, the superscripts refer to the time step 
and the subscripts refer to the computational cell (see Fig. 5). 


Note than in Eq. (10), the change of momentum with 
respect to time is given by the difference of cell edge fluxes. In 
the second equation, neither the velocity nor the pressure are 
written in a conservative form, because the cell-edge difference 
fluxes, that is - p.^) and - (F^)y_i> are 

multiplied by (1/pp and ((F^jp, respectively. 

In particular, for the velocity term. Fig. 5 and Eq. (11) 
indicate that the flux out of the right side of cell / and which 
enters cell /-t 1, is However, the flux into 

cell 0' -t- 1), from cell/, is ((^x)/+i(^x)/-n/ 2 )- Therefore, at 
the interface between two cells, the outgoing flux is different 
than the incoming one and it is for this reason that the 
numerical scheme in Eq. (1 1) is called “non-conservative.” 

It has been shown that the finite difference form of the 
Navier-Stokes equations is frequently unstable, even when the 
equations are written in conservative forms (Refs. 4, 36-38). 
However, these equations can be made stable, by replacing the 
central differences scheme by the “upwind” or “donor-cell” 
differencing scheme (Ref. 39). In this method, the boundary 
values for any dependent variable are a function of the 
direction of the flow velocity. For example, labeling this 
dependent variable by Q, the upwind (donor-cell) scheme 
gives: 

g; (^x);.i/2>o (12) 


ePi if (i;)Pi/2<0 (13) 

This is illustrated in Figs. 6a and 6b. 



V. Simplification of the Governing 
Equations for Wind Modeling 

The Navier-Stokes equation, as written in Section 111, is 
very complicated and its numerical solution would require a 
great deal of computer time and money. Therefore, in 
constructing a numerical scheme, it is useful to see if the 
equations can be simplified in any way. For example, in the 
problem of wind prediction, by using order-of-magnitude 
analysis one can use wind speed and direction measurements 
to check if any of the force and/or acceleration terms are 
significantly smaller than the other terms. As indicated in 
Section I, such measurements do exist for the Goldstone area 
and this data can be used to answer the following questions: 
(1) Is it necessary to keep the temporal acceleration (transient 
term) in the Navier-Stokes equation? (2) Are all the convective 
acceleration terms significant, or can the flow in one direction 
be neglected in comparison with the other two? (3) Can 
viscous effects be neglected and the problem be treated as 
inviscid flow? (4) How can thermal heating and stratification 
be entered into the equations? (5) Are the flow conditions 
laminar or turbulent? This section attempts to answer these 
questions, based on the Goldstone data and on other 
measurements. 

A. Variations in Wind Speed and Direction 
at Goidstone 

The data taken at Goldstone were processed by calculating 
hourly averages of wind speed and direction. Changes in this 
average are illustrated in Fig. 7 for an arbitrarily selected week 
in October 1974. A trend of stronger afternoon winds may be 
discerned, but this trend is not reproducible enough to be 
incorporated in any numerical model. It seems that on a 
diurnal basis there is no steady state for wind velocity, and the 
best approach is to accept an hourly quasi steady state, with 
values given by the averages processed from the Goldstone 
data. Although this means that such a model will have to be 
run for every hour, the model itself will be much simpler to 
run than a transient model and thus the cost per run will be 
lower. The usefulness of this approach has been recognized by 
other investigators who have built models based partially on 
this assumption (Ref. 40). An additional advantage of this 
quasi-steady state approach is that temporary variations in 
turbulent eddy structures will be averaged out; this will also 
contribute in decreasing the complexity of the model. (A note 
on turbulence and its effects appears in Section C.) 

B. Convective Accelerations in Wind Fields 

In the rectangular coordinate system, there are three 
mutually perpendicular velocity components: F^, F^, and V^. 
Spatial changes in these components (e.g., 9F^/9x, dVy/dx, 
dV^ldy) multiplied by the velocities themselves make up the 
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convective accelerations. The Goldstone data indicate that 
wind magnitudes do change from location to location, but that 
velocity changes with height area are, in general, much smaller 
than those in the horizontal plane (see Fig. 8). Such an 
observation is also supported by measurements taken at the 
Savannah River Plant, Georgia (see Fig. 9)(Ref. 41). Note, 
however, that this data represents speed, a scalar, and not 
velocity, a vector. Indeed, Fig. 10 indicates that the wind 
direction changes with height, this change in direction remain- 
ing nearly constant with time. Other data, however, indicate a 
logarithmic distribution of the horizontal velocity components 
(Refs. 8, 42, 43 and 57). 


— U* 



(14) 


In this formula, is the mean speed (m/s) at height Z 
above the ground, U* is the friction velocity^ (m/s), K is 
von Karman’s constant (—0.04), and is the roughness 
length (m). This last quantity is terrain dependent and is 
related to the height of the roughness elements in a given 
terrain, such as trees, buildings, crops, etc. Examples of 
are: 


Terrain 

Roughness length 

Zo,m 

Sand, desert 

0.0003 

Grass 

0.003 - 0.01 

Agricultural crops 

0.04 - 0.20 

Forests 

1.0 -6.0 

City 

1.5 

Snow 

0.00005 - 0001 


It is possible that these conflicting data sets result from 
different wind velocity profiles occurring at different wind 
speeds, and that these changes may also depend on the surface 
roughness. In the model proposed, the acceleration normal to 
the terrain will be taken into account so as to allow the model 
to exhibit terrain sheltering and some aspects of flow 
separation (Ref. 40). 


boundary layers is the problem of turbulence, one of the most 
complex in fluid dynamics. The modeling of turbulence 
in any numerical code is not only difficult, but of uncertain 
accuracy also. Because of this, at least in the preliminary 
program, turbulence effects will not be taken into account. 

D. Thermal Heating 

It was mentioned in Section II that body forces are 
Important in the setup of the Navier-Stokes equations. In the 
problem of air flow, body forces due to gravitational effects 
must be coupled with buoyancy effects. The latter are due to 
thermal heating or, more precisely, to differences in the 
temperature between the air close to the ground (warmer air) 
and air higher up (cooler air). These differences in temperature 
cause differences in density and thus buoyancy forces. The 
Boussinesq parameter accounts for these forces and also 
incorporates the effect of gravity. This parameter can be 
written in either of the two forms: 

ij g = -p(r-r)g (15) 

In these expressions, g is the gravitational acceleration, 
and are reference density and temperature, respectively, 
and d is the volumetric expansion coefficient. 

Note that the second form involves temperatures, which 
must be considered as unknowns for the momentum and 
continuity system of equations; this implies that an additional 
equation, the energy equation, must be added to the system. 

It is not immediately apparent how important the effect of 
density stratification is on the prediction of wind speeds. Most 
wind model studies, (e.g.. Ref. 44) do not compare predicted 
results that include stratification with results that do not. As a 
compromise between accuracy and computational complexity, 
the effect of density stratification will be retained, but instead 
of adding the energy equation to our previous system, the new 
model proposed will incorporate a subroutine giving the 
explicit dependence between density and temperature. The 
advantages of such an approach are also discussed in Ref. 45. 


C. Viscous Effects, Turbulence 

The problems of viscous effects are tied to the velocity 
profiles. If normal accelerations are to be allowed as indicated 
in the previous paragraph, then a full three-dimensional model 
is needed. Such a model must take into effect the influence of 
viscosity, i.e., the boundary layer that is developed between 
the main air flow and the ground. Tied to the nature of 

2 

U* is a measure of surface shear stress. 


E. Simplified System of Equations 

As a result of the above discussions, the system to be solved 
is composed of the following equations: 

Continuity: 


0F dV dV 
dy dz 


= 0 


(16) 


158 



Momentum: 


dV 

7 

* 0a: 


+ V 


0F 0F 

— - + V — - 
by z 02 


P 0a: p 




tions have been worked out by a group of investigators led by 
D. B. Spalding (Refs. 34, 47 and 48). 

The main variables in these procedures are the velocities 
and the pressures. The domain of interest is overlaid by a 
rectangular mesh which can contain cells of different size. The 
differential equations are then approximated by finite differ- 

(17) 

ences applied to this mesh. Figure 1 1 illustrates the placement 
of the variables within the computation grid. 


bV.. 


bV, 






^ bx y by ^ dz 


P P\ dx^ 0y2 9^2 i 


bV bV bV 

V — -+ V — -+ V — - 
X bx y by z 9z 


Note that in this scheme, the variables are staggered, mean- 
ing that the pressure and velocities are stored at different 
locations; the pressure is stored at the intersection of the grid 
lines, e.g., at P, and the velocities at the locations indicated by 
the arrows, i.e., midway between the grid lines. (In this figure, 
(18) the velocity components for point P are and Vy). 

The reason for using a “staggered-grid” arrangement is that 
momentum fluxes can be easily computed from the integra- 
tion of the velocity components over the areas .<4^ andyl^ , and 
such integrations form the basis of these methods (Ref. 62). 


+i|2+it 

P 0Z P 


0^7 0^F 0^f\ 

- + -+ ^ I 

0X^ 0/ 0Z^ / 



(19) 


This system contains four equations and four unknowns: V^, 
Vy, V^, p. As it stands, the system also needs explicit 
relationships between the density and temperature: 


From the integration of the momentum differential equa- 
tions, and by using a combination of central and upwind 
difference schemes, algebraic equations are obtained involving 
the velocity components and/or any scalar quantity that might 
be needed such as effluent concentration, radiation flux, etc. 
Together with the finite difference form of the continuity 
relation, this set of equations is solved by a semi-imphcit^ 
method. 


P = piT) (20) 

and the viscosity and temperature 

p=p(7) (21) 

A final point worth mentioning is the problem of compres- 
sibility. In the equations presented above, air is treated as an 
incompressible fluid. Airflows at moderate temperatures, say 
between 50 and 100°F, and at speeds less than about 100 m/s 
(225 mph), can be considered incompressible (Ref. 46). Since 
these limits are well beyond the situations expected in this 
study of winds, the numerical analysis will be based on an 
incompressible flow. 


VI. Steady-State Models 

The differential equations given in the previous section 
describe steady, fully developed, three-dimensional flows. 
Effective numerical procedures for the solution of these equa- 


The steps in this method are first to estabUsh guessed 
pressures (either as pure guesses or as values from a previous 
calculation sweep), and then to obtain a field of intermediate 
velocities. In general, these intermediate velocities do not 
satisfy the continuity equation and a pressure correction term 
must then be introduced. The process of finding this pressure 
correction field reduces to solving a Poisson equation, a prob- 
lem for which there exist many rapid and economical methods 
of solution (Ref. 49). 

The algorithm used in these papers has been apphed to 
problems such as that of laminar flow over a slab-sided “build- 


^The terms explicit and implicit appear quite often in the description 
of solutions for partial differential equations. Explicit methods are 
those in which the unknown quantity, say a velocity at a time t, V, is 
determined in terms of previously found quantities, say V^~^, 
V^~^, .... In the implicit methods, fivo or more values of the 

unknown quantity are specified in terms of known values on a 
previous data line. For example, the values of say the velocity at time 
t and at different locations can be described as a function of the 
velocities at these locations at a time (r-1). Implicit methods are often 
used in time-dependent heat conduction problems. 
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ing” (Ref. 48), flow with heat transfer in three-dimensional 
ducts (Ref. 50), and boundary layer flows (Ref. 5). The results 
are quite good and the methods appear economical, although 
no figures are given. 

There exist some restrictions imposed on these steady pro- 
grams. These limitations are of three kinds: 

(1) All the papers cited above deal with problems involving 
low Reynolds number laminar flows. At higher Rey- 
nolds numbers, which might be expected in wind fields, 
turbulence may become a factor. In this case, one 
would have to incorporate a turbulence model such as 
those outlined by Gawain and Pritchett (Ref. 52) or by 
Deardorff (Ref. 53). 

(2) The boundaries of the domain of integration are rigid 
and impermeable to fluid flow. This simplifies the 
boundary conditions. In a real situation, however, 
topographical variations will affect the wind field at the 
boundaries. 

(3) Finally, all the problems under consideration have dealt 
with regular geometries: tubes, rectangular blocks, etc. 
In the wind flow problem, the topography is highly 
irregular and storing the necessary terrain description 
can be a significant problem. 

None of the difficulties described above represent insur- 
mountable obstacles. As mentioned above, numerical schemes 
that describe turbulence aheady exist in the literature. The 
problem of irregular topography has been attacked by Viecelli 
(Ref. 54) and also by Mason and Sykes (Ref. 55). Extrapola- 
tion formulas can be used to describe the boundary flows; an 
example of this approach appears in Clark (Ref. 56). It seems 
possible, therefore, to build on these works and write an 
accurate wind-field prediction model. 


VII. Conclusions and Plans for Future Work 

The present report has surveyed various numerical methods 
which can be used to predict the wind field over an arbitrary 
terrain, together with an overview of modeling activities by 
various centers. Several conclusions can be drawn from this 
review. 

Writing such a numerical code is feasible. Numerous papers 
have been cited which successfully predict the dispersal of 
pollutants and the velocity vectors over slabUke, cuboidal and 
bell-shaped obstacles. There is no intrinsic and impenetrable 
difficulty in predicting the velocities over an irregular, three- 
dimensional terrain such as one would encounter in real life. 
No claim is made that the problem is simple and that all one 
has to do is to slap together previously worked out techniques; 
however, writing such a numerical code is perfectly feasible. 

Due to the nature of the problem which requires the 
manipulation of large amounts of data, any wind prediction 
numerical model must be economical to run. This economy 
can be gained by various means — the most obvious one being 
good programming. Beyond this purely technical aspect, how- 
ever, the physical data of the problem can also be used to 
reduce the computational load. For example, by assuming an 
incompressible and steady-state flow, the governing equations 
can be simplified considerably and the computational com- 
plexity is reduced. The only limit in this simpUfication process 
is that the resulting physical model must give a wind field 
prediction reasonably close to terrain measurements. 

A numerical wind model along the lines described above is 
in the process of being written in the DSN Advanced Engineer- 
ing and Energy Conservation Group. This model will take 
advantage of the existing Goldstone Wind Data Base. The 
results of the computer program will be checked both against 
terrain measurements and the wind tunnel simulations. 
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Fig. 3a. Two-dimensional case 
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Fig. 2. Computational mesh for one-dimensional velocity derivative 
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Fig. 3b. Three-dimensional case 













Fig. 7. Wind speed vs time (Goldstone data) 


167 









SPEED, fl/sec 



2 3 

SPEED, m/$ 

9. Wind profile at Savannah River Plant, Ga. 


169 


WIND DIRECTION, deg 



HOUR OF THE DAY 


Fig. 10. Comparison of wind directions at different heights 



Fig. 11. The staggered grid used in the steady-state model 
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